#!/usr/bin/env julia
"""
Basic 2D Navier-Stokes Example

This example demonstrates how to use the NSEMSolver.jl package for solving
2D incompressible Navier-Stokes equations with manufactured solutions.
"""

using NSEMSolver
using Printf

function main()
    println("🌊 NSEMSolver.jl - Basic 2D Example")
    println("=" ^ 50)
    
    # Create NSOptions for the simulation
    options = NSOptions(
        N = 4,                    # Polynomial order
        n_block = 3,             # Domain blocks per direction
        nu = 0.01,               # Kinematic viscosity
        tfinal = 1.0,            # Final time
        cfl = 0.5,               # CFL number
        tol = 1e-6,              # Convergence tolerance
        solver = :julia,         # Linear solver backend
        method = :sem,           # Spectral element method
        adaptive_refinement = true,  # Enable adaptive refinement
        refinement_levels = 3,   # Number of refinement levels
        verbose = true,          # Verbose output
        save_history = true      # Save convergence history
    )
    
    println("Configuration:")
    println("  Polynomial order: $(options.N)")
    println("  Domain blocks: $(options.n_block)")
    println("  Viscosity: $(options.nu)")
    println("  Final time: $(options.tfinal)")
    println("  CFL number: $(options.cfl)")
    println("  Solver: $(options.solver)")
    println()
    
    # Solve the Navier-Stokes equations
    println("Starting simulation...")
    start_time = time()
    
    result = solve_navier_stokes_2d(options)
    
    end_time = time()
    
    # Display results
    println()
    println("📊 Simulation Results:")
    println("  Converged: $(result.converged)")
    println("  Iterations: $(result.iterations)")
    println("  Final residual: $(@sprintf("%.2e", result.residual_norm))")
    println("  Total solve time: $(@sprintf("%.3f", result.solve_time)) seconds")
    
    if result.converged
        println("  ✅ Solution converged successfully!")
    else
        println("  ⚠️  Solution did not converge within tolerance")
    end
    
    # Solution statistics
    if result.u !== nothing
        max_u = maximum(abs.(result.u))
        max_v = maximum(abs.(result.v))
        println("  Maximum |u|: $(@sprintf("%.3f", max_u))")
        println("  Maximum |v|: $(@sprintf("%.3f", max_v))")
    end
    
    # Convergence history
    if length(result.convergence_history) > 0
        println("\n📈 Convergence History (last 5 steps):")
        n_show = min(5, length(result.convergence_history))
        for i in (length(result.convergence_history) - n_show + 1):length(result.convergence_history)
            step_num = length(result.convergence_history) - n_show + (i - length(result.convergence_history) + n_show)
            println("    Step $(step_num): $(@sprintf("%.2e", result.convergence_history[i]))")
        end
    end
    
    # Demonstrate basic visualization capabilities
    println("\n🎨 Creating basic visualizations...")
    
    # Plot the complete solution
    println("  📊 Plotting complete solution...")
    plot_ns_solution(result)
    
    # Create individual velocity component plots
    println("  🌀 Plotting velocity components...")
    plot_velocity_field(result, :magnitude)
    plot_velocity_field(result, :u)
    plot_velocity_field(result, :v)
    
    # Plot pressure contours
    println("  💫 Plotting pressure contours...")
    plot_pressure_contours(result, levels=10)
    
    # Plot convergence history if available
    if !isempty(result.convergence_history)
        println("  📈 Plotting convergence history...")
        plot_convergence_history(result, log_scale=true)
    end
    
    # Create visualization summary
    println("  📋 Creating visualization summary...")
    summary = create_visualization_summary(result)
    
    println("  📝 Visualization Summary:")
    println("    Grid size: $(summary["grid_size"])")
    println("    Velocity magnitude range: [$(summary["velocity_magnitude_stats"]["min"]), $(summary["velocity_magnitude_stats"]["max"])]")
    println("    Pressure range: [$(summary["pressure_stats"]["p_min"]), $(summary["pressure_stats"]["p_max"])]")
    
    # Export results for external analysis
    println("  💾 Exporting results...")
    export_vtk(result, "basic_2d_example_result")
    
    println("\n✅ Visualizations created successfully!")
    println("📁 Files generated (if visualization backends available):")
    println("   - Solution plots (PNG format)")
    println("   - basic_2d_example_result.vtu (VTK format for ParaView)")
    
    println("\n🎉 Example completed successfully!")
    return result
end

# Run the example
if abspath(PROGRAM_FILE) == @__FILE__
    main()
end